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Abstract 

A  new  algorithm  is  developed  based  on  the  complete  two-fluid  plasma  model. 
The  model  captures  physics  that  is  not  contained  in  the  more  common  magne¬ 
tohydrodynamic  model.  Recent  progress  on  the  algorithm  development  for  the 
two-fluid  plasma  model  includes  analysis  of  the  simulation  results,  extension  to 
two  dimensions,  and  preliminary  investigation  of  a  higher  order  method.  The 
algorithm  is  first  developed  in  one  dimension  and  applied  to  the  electromag¬ 
netic  plasma  shock.  Results  are  analyzed  to  provide  physical  support  for  the 
numerically  observed  fast  waves.  The  two-fluid  plasma  algorithm  is  extended 
to  multidimensions  including  a  complete  electromagnetic  field  description.  The 
divergence  constraints  on  the  fields  are  satisfied  by  transforming  the  constraints 
into  hyperbolic  equations.  A  preliminary  higher-order  discontinuous  Galerkin, 
Runge-Kutta  method  is  investigated  to  more  accurately  balance  equilibrium 
sources. 


1  Executive  Summary 

This  project  represents  a  three  year  effort  to  develop  a  new  algorithm  for  plasma 
simulations  based  on  the  two-fluid  plasma  model.  Current  plasma  simulation 
algorithms  capable  of  complex  geometries  are  based  on  the  MHD  (magneto¬ 
hydrodynamic)  model.  The  derivation  of  the  MHD  model  involves  several 
assumptions  that  severely  limit  its  applicability.  The  two-fluid  model  only 
assumes  local  thermodynamic  equilibrium  and  is,  therefore,  more  physically 
accurate  and  capable  than  MHD  models.  The  two-fluid  model  is  formulated 
in  conservation  form.  An  approximate  Riemann  solver  is  developed  for  the 
two-fluid  plasma  model  to  compute  the  fluxes  in  a  stable  and  accurate  man¬ 
ner.  Several  methods  are  investigated  to  solve  the  electromagnetic  field  model, 
which  includes  the  source  terms  and  divergence  constraints.  These  methods 
include  a  characteristic-based  algorithm,  a  perfectly  hyperbolic  modification, 
and  a  mixed  potential  formulation.  The  two  plasma  fluids  and  the  electromag¬ 
netic  fields  communicate  through  the  source  terms.  The  fluid  momentum  and 
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energy  equations  have  source  terms  that  depend  on  E  and  B.  The  electromag¬ 
netic  equations  have  source  terms  that  depend  on  V;  and  ve  (Ampere’s  law) 
and  rii  and  ne  (Poisson’s  equation).  Accurately  coupling  the  source  terms  is 
important  both  for  numerical  stability  and  for  modeling  plasmas  where  large 
equilibrium  forces  exist. 

An  algorithm  is  developed  for  the  complete  two-fluid  plasma  model  ini¬ 
tially  in  one  dimension.  The  algorithm  uses  a  Roe-type  approximate  Riemann 
solver  [1]  to  discretize  the  hyperbolic  fluxes  of  the  fluid  model  and  an  upwind 
characteristic-based  solver  for  the  electromagnetic  fields.  Simulations  from  the 
resulting  finite  volume  algorithm  are  benchmarked  against  known  analytical 
results.  Furthermore,  the  algorithm  is  applied  to  the  electromagnetic  plasma 
shock  problem  to  reveal  the  transition  from  gas  dynamic  shock  to  MHD  shock. 
The  results  are  analyzed  to  reveal  the  fast  plasma  waves  that  are  captured  in 
the  two-fluid  plasma  model.  [2] 

The  algorithm  is  extended  to  multiple  dimensions.  The  divergence  con¬ 
straints  can  be  difficult  to  satisfy  with  the  presence  of  current  and  charge 
sources  on  an  arbitrary  computational  grid.  The  divergence  constraints  are 
satisfied  by  reformulating  Maxwell’s  equations  to  include  correction  potentials. 
The  approach  involves  coupling  the  divergence  constraint  equations  with  the 
time-dependent  field  equations  to  form  a  perfectly  hyperbolic  equation  set.  [3] 
The  finite  volume  algorithm  is  applied  to  a  two-dimensional  plasma  shock  prob¬ 
lem  which  is  similar  to  a  field  reversed  configuration  (FRC)  being  investigated 
at  the  Air  Force  Research  Laboratories  (FRC  magnetic  implosions).  The  sim¬ 
ulations  demonstrate  the  importance  of  accurately  modeling  the  source  terms. 
The  FRC  simulation  must  balance  large  magnetic  field  forces  with  large  plasma 
pressure  forces.  A  preliminary  algorithm  is  developed  using  a  discontinuous 
Galerkin  method.  [4]. 

This  project  resulted  in  the  Master  of  Science  theses  for  Chris  Aberle  and 
John  Loverich  titled  “Algorithm  for  Solving  Colocated  Electromagnetic  Fields 
with  Sources”  and  “A  Finite  Volume  Algorithm  for  the  Two-fluid  Plasma  Sys¬ 
tem  in  One  Dimension”,  respectively.  Archival  journal  and  conference  papers 
were  published  reporting  on  the  work  from  this  project.  Archival  journal  papers 
are  U.  Shumlak,  J.  Loverich,  “Approximate  Riemann  solver  for  the  two-fluid 
plasma  model”,  Journal  of  Computational  Physics  187,  620-638  (2003)  and  “A 
discontinuous  Galerkin  method  for  the  full  two- fluid  plasma  model”,  Computer 
Physics  Communications  169  251-255  (2005).  This  project  was  performed  by 
Prof.  Uri  Shumlak  and  graduate  students  Chris  Aberle,  Ammar  Hakim,  and 
John  Loverich. 


2  Project  Description 

Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force, 
some  of  which  have  dual-use  potential.  These  applications  include  portable 
pulsed  power  systems,  high  power  microwave  devices,  drag  reduction  for  hyper- 


sonic  vehicles,  advanced  plasma  thrusters  for  space  propulsion,  nuclear  weapons 
effects  simulations,  radiation  production  for  counter  proliferation,  and  fusion 
for  power  generation.  Several  of  these  applications  are  specifically  mentioned 
in  the  New  World  Vistas  Report  from  the  USAF  Scientific  Advisory  Board.  [5] 
In  general,  plasmas  fall  into  a  density  regime  where  they  exhibit  both  collective 
(fluid)  behavior  and  individual  (particle)  behavior.  The  intermediate  regime 
complicates  the  computational  modeling  of  plasmas. 


2.1  Kinetic,  PIC,  and  MHD  Models 


Plasmas  may  be  most  accurately  modeled  using  kinetic  theory.  The  plasma  is 
described  by  distribution  functions  in  physical  space,  velocity  space,  and  time, 
/(x,  v,t).  The  evolution  of  the  plasma  is  then  modeled  by  the  Boltzmann 
equation. 


dfa 

dt 


+  VQ 


|  +  -(E  +  vaxB) 
ox  ma 


dfa 

dv 


dfa 


dt 


collisions 


(i) 


for  each  plasma  species  a  =  ions,  electrons.  The  Boltzmann  equation  cou¬ 
pled  with  Maxwell’s  equations  for  electromagnetic  fields  completely  describe 
the  plasma  dynamics.  [6-8]  However,  the  Boltzmann  equation  is  seven  dimen¬ 
sional.  As  a  consequence  of  the  large  dimensionality  plasma  simulations  using 
the  Boltzmann  equation  are  only  used  in  very  limited  applications  with  narrow 
distributions,  small  spatial  extent,  and  short  time  durations.  [9, 10]  The  seven 
dimensional  space  is  further  exacerbated  by  the  high  velocity  space  that  is  un¬ 
used  except  for  tail  of  the  distribution  or  energetic  beams.  Boundary  conditions 
are  difficult  to  implement  in  kinetic  simulations. 

Particle  in  cell  (PIC)  plasma  model  apply  the  Boltzmann  equation  to  repre¬ 
sentative  superparticles  which  are  far  fewer  (107)  than  the  number  of  particles 
in  the  actual  plasma  (1020).  [11]  PIC  simulations  have  similar  limitations  as 
simulations  using  kinetic  theory  due  to  statistical  errors  caused  by  the  fewer 
superparticles.  Boundary  conditions  are  also  difficult  to  implement  in  PIC 
simulations. 

The  other  end  of  the  spectrum  in  plasma  model  involves  taking  moments 
of  the  Boltzmann  equation  and  averaging  over  velocity  space  for  each  species 
which  implicitly  assumes  local  thermodynamic  equilibrium.  The  resulting  equa¬ 
tions  comprise  the  two-fluid  plasma  model.  The  two-fluid  equations  are  then 
combined  to  form  the  MHD  model.  [12]  However,  in  the  process  several  ap¬ 
proximations  are  made  which  limit  the  applicability  of  the  MHD  model  to  low 
frequency  and  ignores  the  electron  mass  and  finite  Larmor  radius  effects. 

The  MHD  model  treats  the  plasma  like  a  conducting  fluid  and  assigning 
macroscopic  parameters  to  describe  its  particle-like  interactions.  Plasma  simu¬ 
lation  algorithms  based  on  the  MHD  model  have  been  very  successful  in  model¬ 
ing  plasma  dynamics  and  other  phenomena.  Codes  such  as  MACH2  are  based 
on  arbitrary  Lagrangian/Eulerian  formulations.  [13]  ALE  codes  are  well  suited 
for  simulating  plasma  phenomena  involving  moving  interfaces.  [14]  However, 
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Figure  1:  Dispersion  relation  for  the  fast  and  slow  modes  with  uict  =  1.0.  The  dashed 
lines  represent  the  ideal  MHD  modes. 

ALE  codes  cannot  be  formulated  as  conservation  laws  and  lack  many  of  the 
inherent  conservative  properties.  The  MHD  model  has  been  successfully  im¬ 
plemented  in  conservative  form  to  simulate  realistic  three-dimensional  geome¬ 
tries.  [15, 16] 

A  severe  limitation  of  the  MHD  model  is  the  treatment  the  Hall  effect 
and  diamagnetic  terms.  These  terms  represent  the  separate  motions  of  the 
ions  and  electrons.  The  Hall  effect  and  diamagnetic  terms  also  account  for 
ion  current  and  the  finite  ion  Larmor  radius.  These  effects  are  important  in 
many  applications  such  as  electric  space  propulsion  thrusters:  Hall  thrusters, 
magnetoplasmadynamic  (MPD)  thrusters,  Lorentz  force  thrusters.  The  Hall 
term  is  also  believed  to  be  important  to  electrode  effects  such  as  anode  and 
cathode  fall  which  greatly  affect  many  directly  coupled  plasma  applications. 
Furthermore,  the  Hall  and  diamagnetic  effects  may  be  important  for  hypersonic 
flow  applications.  [17] 

The  Hall  terms  can  be  difficult  to  stabilize  because  they  lead  to  the  whistler 
wave  branch  of  the  dispersion  relation.  The  phase  and  group  velocities  of  the 
whistler  wave  increase  with  frequency.  The  velocities  become  large  even  for 
modest  values  of  the  Hall  parameter.  See  Fig.  1  for  the  dispersion  diagram  for 
a  modest  value  of  the  Hall  parameter  (a ;ct  =  1). 

A  semi-implicit  technique  has  been  applied  to  treat  the  Hall  term.  [18, 19] 

After  the  hyperbolic  terms  of  the  MHD  equations  are  advanced,  the  Hall  terms 
are  treated  independently.  The  conserved  variables  are  then  corrected.  The 
procedure  can  be  computationally  intensive.  The  operator  stencil  uses  5  points 
in  the  sweep  direction  and  3  points  in  each  orthogonal  direction.  The  complete 
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operator  stencil  is  45  points.  The  semi-implicit  method  works  adequately  for 
small  Hall  parameters,  but  becomes  unstable  or  slow  to  converge  for  the  large 
Hall  parameters  often  seen  in  applications. 


2.2  Project  Results 

The  next  step  is  to  extend  the  realistic  geometry  capabilities  of  the  MHD  model 
to  the  more  physically  accurate  two-fluid  model.  The  complexity  of  the  two- 
fluid  model  is  greater  the  MHD  model  but  significantly  less  than  the  kinetic 
model. 

In  this  project  a  new  algorithm  is  developed  that  solves  the  two-fluid  plasma 
model  using  an  approximate  Riemann  solver.  [2]  The  method  tracks  the  wave 
propagation  across  the  domain  based  on  conservation  laws. 

The  two-fluid  plasma  model  captures  the  separate  motion  of  the  electrons 
and  ions  without  the  added  complexity  of  the  kinetic  model.  The  two-fluid 
model  is  derived  by  taking  moments  of  the  Boltzmann  equation  for  each  species. 
The  process  of  taking  moments  eliminates  the  velocity  space  and  yields  rep¬ 
resentative  fluid  variables  (density,  momentum,  energy)  for  each  species.  The 
only  approximation  made  is  local  thermodynamic  equilibrium  within  each  fluid 
and  is,  therefore,  a  generalization  of  the  MHD  model.  The  fundamental  vari¬ 
ables  are  generated  by  taking  moments  of  the  distribution  function. 

The  evolution  of  the  particle  density  of  the  ions  and  electrons  is  expressed  by 
continuity  equations.  The  equations  are  the  zeroth  moment  of  the  Boltzmann 
equation. 

^-v(i)=0  <3) 

where  n*,  ne  are  the  ion  and  electron  number  densities  and  the  particle  fluxes 
are  defined  by  the  partial  current  densities  jj  =  gn,Vj  and  je  =  —eneve  in  terms 
of  the  charges  q,  e  and  fluid  velocities  Vj,  ve  for  the  ion  and  electron  fluids. 

The  first  moment  of  the  Boltzmann  equation  yields  momentum  equations 
for  each  species.  The  momentum  equations  are  written  in  divergence  form. 
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where  E  and  B  are  the  electric  and  magnetic  fields,  pi  and  pe  are  the  ion  and 
electron  partial  pressures,  and  Re;  is  the  electron  to  ion  momentum  transfer 
vector. 


The  second  moment  of  the  Boltzmann  equation  yields  energy  equations  for 
each  species  which  are  expressed  in  divergence  form  for  the  total  energy. 
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where  the  total  energy  is  defined  by 
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where  7  is  the  ratio  of  specific  heats  and  T{,  Te  are  the  ion  and  electron  tem¬ 
peratures.  An  adiabatic  equation  of  state  is  assumed.  The  temperatures  are 
related  to  the  partial  pressures  by  pa  =  naTa  for  a  =  {i,  e}. 

The  equations  that  govern  the  ion  and  electron  fluids  are  rewritten  in  com¬ 
pact,  conservative  form. 

^  +  V-F  =  S  (10) 

where  Q  is  the  vector  of  conserved  fluid  variables,  F  is  the  tensor  of  hyperbolic 
fluxes  ( Fx  +  Gy  +  Hz),  and  S  is  the  vector  containing  the  source  terms.  The 
vector  of  conserved  variables  is 


Q  [rq ,  ne,  jix ,  jiy ,  jiz ,  j£x ,  jey ,  jez ,  S{,  £e]  .  (11) 

for  the  number  densities,  electrical  current  densities,  and  energy  densities.  The 
flux  Jacobian  dF/dQ  for  the  two  fluid  equations  is  constructed  in  the  usual 
way.  The  characteristic  velocities  are  calculated  to  construct  the  approximate 
Riemann  fluxes. 

The  eigenvalues  of  the  flux  Jacobian  give  the  characteristic  velocities. 

^  J**  -J-  jex  jex 

|  en-i  ’  erii  \  3  raj  ’  ene  ’  ene 

The  eigenvalues  for  the  two  fluid  plasma  model  represent  the  combination  of 
the  drift  speeds  and  thermal  speeds  for  the  electrons  and  ions. 

The  electromagnetic  fields  influence  the  motion  of  the  plasma  fluid  through 
the  Lorentz  force  which  is  contained  in  eqns(4)  and  (5).  The  motion  of  the 
plasma  influences  the  evolution  of  the  electromagnetic  fields  through  the  re¬ 
distribution  of  charge  density  and  current  density.  Maxwell’s  equation  govern 
the  evolution  of  the  electromagnetic  fields.  The  charge  density  qrii  —  ene  and 


current  density  (j  =  j ;  +  je)  are  calculated  directly  from  the  two-fluid  equations 
which  couple  the  electromagnetic  fields. 

#=-ViE  <13> 

eo~Q^  =  V  x  B//i°  -  j  (14) 

eoV  ■  E  =  qrii  -  ene  (15) 

V  •  B  =  0  (16) 

The  two-fluid  plasma  model  (including  the  electromagnetic  equations)  can 
also  be  expressed  in  divergence  form. 

^  +  V-F  =  S  (17) 

where  Q  is  the  vector  of  conserved  fluid  variables,  F  is  the  tensor  of  hyperbolic 
fluxes  ( Fx  +  Gy  +  Hz),  and  S  is  the  vector  containing  the  source  terms. 

The  hyperbolic  fluxes  are  discretized  using  a  Roe-type  approximate  Rie- 
mann  solver.  [1]  In  this  method  the  overall  solution  is  built  upon  the  solutions 
to  the  Riemann  problem  defined  by  the  discontinuous  jump  in  the  solution 
at  each  cell  interface.  The  numerical  flux  for  a  first-order  accurate  (in  space) 
Roe-type  solver  is  written  in  symmetric  form  as 


Fi+l/2  2 


\  (Fi+ 1  +  Fi)-~^2lk  ( Qi+1  -  Qi)  |Afc|  rk 


where  rk  is  the  kth  right  eigenvector,  Xk  is  the  kth  eigenvalue,  and  lk  is  the 
kth  left  eigenvector,  evaluated  at  the  cell  interface  (i  +  1/2).  The  values  at  the 
cell  interface  are  obtained  by  a  Roe  average  of  the  neighboring  cells.  Higher 
order  accuracy  is  obtained  by  using  a  minmod  flux  limiter  in  the  numerical 
flux  calculation.  The  algorithm  is  second-order  accurate  in  regions  where  the 
solution  is  smooth  and  first-order  accurate  in  the  vicinity  of  large  gradients  in 
the  solution.  The  flux  calculated  as  above  is  normal  to  the  cell  interface  which 
is  the  desired  orientation  for  applying  the  divergence  theorem  in  a  finite  volume 
method. 

The  source  terms  of  eqn(17)  are  large,  in  general,  since  they  contain  the 
electromagnetic  forces  that  act  upon  the  fluids.  The  homogeneous  fluxes  are 
calculated  for  the  electron  and  ion  fluids  with  the  approximate  Riemann  solver 
described  previously.  An  approximate  Riemann  solver  is  also  used  to  calculate 
the  homogeneous  fluxes  for  Maxwell’s  equations.  The  large  source  terms  make 
the  equation  stiff.  The  source  terms  for  all  of  the  equations  are  then  solved  self- 
consistently  with  an  implicit  treatment  to  alleviate  the  stiffness  of  the  problem. 


8 


Figure  2:  Two-fluid  solution  of  the  electromagnetic  plasma  shock  problem  for  four 
values  of  the  Larmor  radius  77,.  The  transition  from  gas  dynamic  shock  to  MHD 
shock  is  evident.  Fast  moving  waves  are  present  for  intermediate  values  of  ri. 


where  5n+1  is  the  source  term  evaluated  at  the  n  + 1  time  step.  Previously,  the 
source  term  is  expanded  in  a  Taylor  series  to  first  order  and  a  Newton  iteration 
is  required  to  solve  the  equation.  This  step  has  been  eliminated  by  using  an 
implicit  time  approach  to  rewrite  the  equation  as 


dQ 

dt 


71+1 

sources 


=  sn+1 


=  sn  + 


dS_dQ 
dQ  dt 


At  = 


This  formulation  does  not  require  a  Newton  iteration. 


(20) 


2.2.1  One-Dimensional  Algorithm 

The  described  algorithm  is  used  to  simulate  the  coplanar  Riemann  problem 
to  test  for  proper  wave  capturing  behavior,  Langmuir  oscillations  to  test  the 
electric  field  component  in  the  Lorentz  force,  and  the  electromagnetic  plasma 
shock. 

The  two-fluid  algorithm  is  benchmarked  to  the  electromagnetic  plasma 
shock  problem  often  used  for  MHD  algorithms  as  in  Ref.  [20].  The  MHD  shock 
problem  is  a  good  test  of  the  two-fluid  algorithm’s  ability  to  model  MHD  like 
problems.  MHD  problems  are  characterized  by  strong  coupling  of  the  two-fluid 
source  terms  which  makes  for  a  stiff  system  of  equations.  In  this  benchmark 
case,  it  is  shown  that  the  algorithm  can  be  used  to  model  MHD  like  problems. 
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Figure  3:  Fourier  analysis  of  the  numerical  results  for  a  small  amplitude  electro¬ 
magnetic  plasma  shock  compared  to  the  analytical  dispersion  relations  for  plasma 
waves  for  77,  —  0.1.  Analytical  dispersion  relations  for  the  L  mode,  R  mode,  and 
Alfven  mode  waves  are  plotted  for  comparison.  The  low  frequency  portion  of  the 
lower  branch  of  the  R  mode  is  the  whistler  wave.  The  agreement  of  the  numerical 
and  analytical  results  supports  the  physical  basis  for  the  fast  waves  in  the  two-fluid 
solutions. 

The  comparison  is  made  between  the  ideal  MHD  solution  and  a  two-fluid  solu¬ 
tion  with  equivalent  initial  conditions.  The  two- fluid  model  allows  the  flexibility 
to  choose  the  desired  Larmor  radius  ri  (zero  in  MHD),  the  ionic  charge  state 
(unity  in  MHD),  and  the  electron  to  ion  mass  ratio  (zero  in  MHD).  For  the 
simulations  presented  here,  the  ionic  charge  state  is  set  to  one,  and  the  electron 
to  ion  mass  ratio  is  set  to  5.45  x  10-3  (hydrogen  plasma).  The  Larmor  radius 
is  varied  while  holding  the  Alfven  and  sound  speeds  constants. 

Simulation  results  are  presented  in  Fig.  2  comparing  the  results  for  four 
values  of  the  Larmor  radius.  As  the  Larmor  radius  decreases,  the  two-fluid 
solution  approaches  the  MHD  solution.  The  algorithm  is  able  to  capture  the 
MHD  shock  in  the  limit  of  77,  — >  0,  the  gas  dynamic  shock  in  the  limit  of 
77  — >  00,  and  two-fluid  (Hall  effect)  physics  for  intermediate  values.  The  MHD 
model  assumes  77,  =  0  where  the  electron  and  ion  fluids  are  tightly  coupled  to 
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the  magnetic  field.  The  plasma  completely  decouples  from  the  magnetic  field 
when  rj,  =  oo.  The  transition  from  gas  dynamic  shock  to  MHD  shock  is  a 
smooth  transition  that  shows  faster  moving  waves  for  intermediate  values  of 
the  Larmor  radius  and  has  not  been  previously  studied. 

Simulation  results  reveal  the  transition  from  gas  dynamic  shock  to  MHD 
shock  is  a  smooth  transition  that  shows  fast  moving  waves  for  intermediate 
values  of  the  Larmor  radius  and  had  not  been  previously  studied.  The  character 
of  these  fast  waves  is  studied  by  decreasing  the  shock  strength  so  nonlinear 
effects  become  small  and  linear  wave  analysis  can  be  applied.  The  simulation 
results  are  Fourier  analyzed  to  form  a  dispersion  plot  which  is  compared  to  the 
analytic  dispersion  relations  from  plasma  wave  theory.  For  the  electromagnetic 
plasma  shock  problem  the  wave  vector  is  perpendicular  to  the  discontinuity 
and  parallel  to  the  longitudinal  magnetic  field.  The  plasma  dispersion  relation 
for  this  case  yields  the  left  and  right  circularly  polarized  waves  (L  mode  and  R 
mode),  in  addition  to  the  slower  Alfven  waves. 


c2k2 


=  1 


u. 


pe 


U) * 


U)  (iv  ±  U)ce) 


(21) 


where  u>pe  =  y/nee2/mee0  is  the  electron  plasma  frequency  and  u>cc  —  eB/me 
is  the  electron  cyclotron  frequency.  The  simulation  and  analytical  results  are 
shown  in  Fig.  3.  The  low  frequency  portion  of  the  lower  branch  of  the  R  mode 
is  the  whistler  wave.  The  agreement  of  the  numerical  and  analytical  results 
supports  the  physical  basis  for  the  fast  waves  in  the  two-fluid  simulations. 


2.2.2  Multi-Dimensional  Algorithm 

The  formulation  of  the  approximate  Riemann  solver  for  the  fluid  equations  al¬ 
lows  for  a  straight  forward  extension  to  multiple  dimensions.  The  extension 
from  one  to  two  dimensions  introduces  the  complexity  of  extending  the  algo¬ 
rithm  to  multiple  dimensions.  The  extension  from  two  to  three  dimensions  only 
requires  additional  compute  resources  but  does  not  add  additional  complexity. 
Therefore,  a  large  effort  is  concentrated  on  the  two  dimensional  algorithm. 

The  divergence  constraints  on  the  magnetic  and  electric  fields  are  satisfied 
automatically  in  one  dimension  and  the  electromagnetic  fields  are  solved  using 
an  upwind  characteristic-based  method.  However,  in  multiple  dimensions  the 
divergence  constraints  can  be  difficult  to  satisfy  with  the  presence  of  current 
and  charge  sources  on  an  arbitrary  computational  grid.  The  divergence  con¬ 
straints  are  satisfied  by  reformulating  Maxwell’s  equations  to  include  correction 
potentials. 

The  approach  involves  coupling  the  divergence  constraint  equations  with  the 
time-dependent  field  equations  to  form  a  perfectly  hyperbolic  equation  set.  [3] 
The  field  equations  are  expressed  as 

SB 

— +  VxE  +  7V^  =  0,  (22) 
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(23) 


(24) 

(25) 


where  <p  and  ip  are  the  electric  and  magnetic  correction  potentials  or,  more 
formally,  Lagrange  multipliers  which  vanish  at  the  domain  boundaries.  The 
method  more  accurately  predicts  the  propagation  speed  of  the  waves;  how¬ 
ever,  the  method  can  overestimate  the  Lorentz  force  caused  by  charge  sepa¬ 
ration.  The  implementation  illustrates  the  importance  of  tightly  coupling  the 
field  solver  to  the  fluid  solver. 


Open  or  absorbing  boundaries  can  be  troublesome  for  finite  volume  meth¬ 
ods  applied  to  electromagnetic  fields  because  the  impedance  mismatch  causes 
spurious  reflections  and  oscillations.  The  problem  has  been  solved  by  imple¬ 
mentation  a  PML  (perfectly  matched  layer)  boundary  condition  [21,22]  which 
is  common  in  PIC  codes. 


The  two  dimensional  algorithm  has  been  applied  to  the  cylindrical  elec¬ 
tromagnetic  plasma  shock  problem  and  the  collisionless  magnetic  reconnec¬ 
tion  problem.  The  simulation  results  for  the  two-dimensional  electromagnetic 
plasma  shock  are  shown  in  Fig.  4.  The  fast  waves  have  already  left  the  do¬ 
main  by  this  time.  The  discontinuous  jump  in  the  out-of-plane  magnetic  field 
is  supported  by  an  azimuthal  current  carried  primarily  by  the  electron  fluid. 
The  electron  fluid  undergoes  a  Kelvin-Helmholtz  instability  which  broadens 
the  current  layer.  (This  behavior  may  be  important  for  shock  degradation  in 
hypersonic  vehicles.)  The  simulation  results  from  the  collisionless  magnetic 
reconnection  problem  are  shown  in  Fig.  5.  The  simulation  demonstrates  an 
important  numerical  effect. 

As  described,  the  source  terms  make  the  equations  stiff  and  must  be  treated 
implicitly  to  relax  the  stiffness.  The  relaxation  works  well  for  transient  prob¬ 
lems.  However,  for  problems  involving  equilibrium  conditions,  the  equilibrium 
force  balance  artificially  diffuses  resulting  in  a  loss  of  the  original  equilibrium. 
The  simulation  ultimately  reduces  to  a  trivial  equilibrium. 

Magnetic  reconnection  is  the  process  where  plasma  currents  coalesce  and 
magnetic  fields  reconnect.  The  process  is  demonstrated  by  a  plasma  sheet  car¬ 
rying  an  initially  uniform  current.  The  current  produces  oppositely  directed 
magnetic  fields  above  and  below  the  plasma  sheet.  The  plasma  current  sup¬ 
ports  the  jump  in  the  magnetic  field,  and  the  configuration  is  in  force  balance 
(equilibrium)  though  large  forces  are  present.  If  the  plasma  has  collisions  (finite 
resistivity),  the  current  sheet  breaks  into  small  filaments  through  a  magnetic 
tearing  instability.  [23]  However,  even  perfectly  conducting  plasmas  undergo  a 
similar  instability  that  is  due  to  two-fluid  effects  instead  of  resistivity.  The  cur¬ 
rent  sheet  again  breaks  into  current  filaments.  The  magnetic  field  reconnects 
“through”  the  remaining  plasma  sheet.  The  process  is  important  in  many  ap¬ 
plications  including  space  physics.  [24]  However,  it  is  a  difficult  problem  to 
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Figure  4:  Two-dimensional  electromagnetic  plasma  shock  simulation  showing  the 
magnitude  of  in-plane  electron  current  density.  The  azimuthal  current  supports  the 
discontinuous  jump  in  the  out-of-plane  magnetic  field.  The  current  layer  is  seen  to 
broaden  due  into  a  turbulent  layer  due  to  a  Kelvin-Helmholtz  instability. 
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Figure  5:  Magnetic  field  contours  for  the  collisionless  magnetic  reconnection  problem 
solved  with  the  finite  volume  method.  Shown  is  the  magnitude  of  the  transverse  mag¬ 
netic  field  Bx  (parallel  to  the  current  sheet).  The  equilibrium  jump  in  the  magnetic 
field  has  diffused  significantly. 


Figure  6:  Magnetic  field  contours  for  the  collisionless  magnetic  reconnection  problem 
solved  with  a  higher-order  method.  Shown  is  the  magnitude  of  the  transverse  mag¬ 
netic  field  Bx  (parallel  to  the  current  sheet).  The  equilibrium  jump  in  the  magnetic 
field  is  maintained. 
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Figure  7:  Comparison  of  the  transverse  magnetic  field  at  x  =  0.5  computed  using  the 
finite  volume  method  and  the  higher  order  method.  The  higher  order  method  better 
preserves  the  equilibrium  jump  in  the  magnetic  field. 
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model  and  provides  a  rigorous  test  for  the  algorithm.  Simulations  of  the  two- 
dimensional  collisionless  reconnection  problem  are  performed  using  the  algo¬ 
rithm  described  previously  with  the  implicit  source  treatment.  Figure  5  shows 
the  simulation  results.  Shown  are  contours  of  the  magnitude  of  the  transverse 
magnetic  field  Bx  (parallel  to  the  current  sheet).  The  equilibrium  jump  in  the 
magnetic  field  has  diffused  significantly. 

Preliminary  results  using  higher-order  methods  have  solved  the  problem.  [4] 
Figure  6  shows  the  same  simulation  of  the  two-dimensional  collisionless  recon¬ 
nection  problem  using  an  algorithm  that  is  second-order  in  space  and  third- 
order  in  time  implemented  for  a  uniform  computational  mesh.  Figure  7  com¬ 
pares  Bx  through  the  domain  midplane  as  computed  by  the  previous  method 
and  the  higher-order  method.  The  higher-order  method  preserves  the  equilib¬ 
rium  jump  in  the  magnetic  field.  The  higher-order  algorithm  uses  a  discontin¬ 
uous  Galerkin  method  with  a  TVD  Runge-Kutta  time  advance.  [25-27] 

The  conserved  variables  of  the  two-fluid  plasma  model  are  modeled  with  a 
set  of  basis  functions,  {u^}.  The  governing  equations,  expressed  as  eqn  (17), 
are  multiplied  by  each  basis  function  and  integrated  over  the  mesh  cell  volume 
Q.  An  integral  equation  is  generated  for  each  basis  function. 

[  vh^dV  +<f  vhF  ■  dS  -  [  F  •  VvhdV  =  f  vhSdV  (26) 
Jo.  dt  Jdu  J  n  Jn 

where  the  divergence  theorem  has  been  applied  to  the  second  term.  The  volume 
and  surface  integrals  are  replaced  with  Gaussian  quadrature.  The  flux  F  is 
computed  with  the  approximate  Riemann  solver  with  a  limiter  applied  directly 
to  the  conserved  variables  to  get  high  resolution.  The  source  terms  are  described 
by  the  basis  functions  and  are,  therefore,  the  same  order  accurate  as  the  solution 
variables.  This  satisfies  the  higher  order  accuracy  requirement  to  preserve  the 
equilibrium  balance  between  the  divergence  of  the  flux  and  the  source. 

For  the  preliminary,  two-dimensional,  second-order  accurate  algorithm,  a 
linear  set  of  basis  functions  are  used. 

W  =  =  (27) 

where  the  center  of  the  mesh  cell  is  located  at  (xij,  yij )  and  extends  Ax  by  Ay. 
The  conserved  variables  Q  are  defined  as 

Q  =  Qo  +  QxV x  +  QyVy  (28) 

within  each  mesh  cell.  Update  equations  for  the  coefficients  for  each  conserved 
variable  are  found  directly  from  eqn  (26)  applied  to  each  mesh  cell. 

The  temporal  evolution  is  determined  with  a  Runge-Kutta  method. 

The  preliminary  discontinuous  Galerkin  method  has  greatly  reduced  the 
diffusion  of  the  equilibrium  forces.  Figure  6  shows  the  same  simulation  of  the 
two-dimensional  collisionless  reconnection  problem  using  an  algorithm  that  is 
second-order  in  space  and  third-order  in  time.  Figure  7  compares  Bx  through 
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the  domain  midplane  as  computed  by  the  previous  method  and  the  higher 
order  method.  The  higher  order  method  preserves  the  equilibrium  jump  in  the 
magnetic  field. 
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